Set up workspace, load data, and load required packages.
rm(list=ls(all=TRUE))
results <- read.csv("Data/Urchin_ClimateChange_Data.csv", header=T, na.strings="NA")
library("lme4") #linear mixed modeling
library("lmtest") #linear mixed modeling
library("lmerTest") #calculate p-values in models
library("car") #levenes tests
library("emmeans") #post-hoc tests
library("lsmeans") #post-hoc tests
library("plotrix") #calculate std.error
library("rstatix") #use for calculating effect size for power test - not installing properly
library("plyr") #splitting, applying, and combining data
library("dplyr")
library("cowplot") #grid plotting and arrange plots
library("effects") #plot effects of modeling
library("MuMIn")
library("car") #levenes tests, Anova, qqPlots,
library("tidyverse")
library("stats")
library("onewaytests") #allows for Welch test with unequal variances (spine lengths)
library("stats")
library("effsize")
library("FSA")
library("rcompanion")
library("broom") #turns model output into data frames
library("knitr")
library("hms")
library("tidyverse")
library("pwr") #run power tests
library("effsize")
1) Build and run model
2) Check for normality of residuals
3) Check for homogeniety of variance of residuals
4) Look at model summary
5) Run anova to check for significance of main effects
EnviroMean <-
results %>%
select("Day", "Temperature", "pH", "Temp","pHspec", "pCO2out", "AlkTotal") %>%
filter(Day >= "1", Day <= "123") %>%
drop_na(Day) %>%
group_by(Temperature, pH) %>%
summarise(mean_T = mean(Temp, na.rm = T),
se_T = sd(Temp, na.rm = T)/sqrt(6),
min_T = min(Temp, na.rm = T),
max_T = max(Temp, na.rm = T),
mean_pH = mean(pHspec, na.rm = T),
se_pH = sd(pHspec, na.rm = T)/sqrt(6),
mean_pCO2 = mean(pCO2out, na.rm = T),
se_pCO2 = sd(pCO2out, na.rm = T)/sqrt(6),
mean_AlkTotal = mean(AlkTotal, na.rm = T),
se_AlkTotal = sd(AlkTotal, na.rm = T)/sqrt(6))
EnviroMean
## # A tibble: 4 × 12
## # Groups: Temperature [2]
## Temperature pH mean_T se_T min_T max_T mean_pH se_pH mean_pCO2 se_pCO2
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ambient acidifi… 24.8 0.887 22.1 27.4 7.64 0.0102 1083. 65.3
## 2 ambient ambient 24.8 0.886 22.1 27.4 7.96 0.0147 481. 22.6
## 3 heated acidifi… 26.7 0.827 23.8 29.4 7.66 0.0129 1107. 77.2
## 4 heated ambient 26.7 0.818 23.8 29.3 7.96 0.0141 524. 20.9
## # … with 2 more variables: mean_AlkTotal <dbl>, se_AlkTotal <dbl>
Environmental <-
results %>%
select("Day", "Temp", "pH", "Temperature","pHspec","Header", "TankID") %>%
filter(Day >= "1", Day != "126") #use only days from start of experiment when desired conditions were reached
#temp between high and ambient temp
tempmod<-aov(Temp ~ Temperature, data=Environmental)
summary(tempmod) #temperature was different between high and ambient treatment
#pH between high and ambient co2
summary(aov(pHspec ~ pH, data=Environmental)) #pH was different between high and ambient co2 treatments
##HEADERS
#temp between headers in high
aov1<- Environmental %>%
filter(Temperature=="heated")
modaov1<-aov(Temp ~ as.factor(Header), data=aov1)
summary(modaov1) #not different between headers
#temp between headers in ambient temp
aov2<- Environmental %>%
filter(Temperature=="ambient")
summary(aov(Temp ~ as.factor(Header), data=aov2)) #not different between headers
#pH between headers in acidified
aov3<- Environmental %>%
filter(pH=="acidified")
summary(aov(pHspec ~ as.factor(Header), data=aov3)) #not different between headers
#pH between headers in ambient pH
aov4<- Environmental %>%
filter(pH=="ambient")
summary(aov(pHspec ~ as.factor(Header), data=aov4)) #not different between headers
##TANKS
#temp between tanks in high temp
aov1<- Environmental %>%
filter(Temperature=="heated")
summary(aov(Temp ~ as.factor(TankID), data=aov1)) #not different between heated tanks
#temp between tanks in ambient temp
aov2<- Environmental %>%
filter(Temperature=="ambient")
summary(aov(Temp ~ as.factor(TankID), data=aov2)) #not different between ambient tanks
#pH between tanks in acidified
aov3<- Environmental %>%
filter(pH=="acidified")
summary(aov(pHspec ~ as.factor(TankID), data=aov3)) #not different between tanks in acidified
#pH between headers in ambient pH
aov4<- Environmental %>%
filter(pH=="ambient")
summary(aov(pHspec ~ as.factor(TankID), data=aov4)) #not different between headers
| Treatment Comparison | Pr(>F) |
|---|---|
| Heated v. Ambient | <0.001 |
| Acidified v. Ambient | <0.001 |
| Ambient temp headers | 0.998 |
| Heated headers | 0.842 |
| Ambient pH headers | 0.129 |
| Acidified headers | 0.150 |
| Ambient temp tanks | 1.000 |
| Heated tanks | 0.999 |
| Ambient pH tanks | 0.675 |
| Acidified tanks | 0.888 |
Check assumptions of treatment conditions model
#check assumptions
# 1. Normality of residuals
qqPlot(residuals(tempmod)) #very much not normal
shapiro.test(residuals(tempmod)) # nope
hist(residuals(tempmod)) #nope...
Check using nonparametric Kruskal-Wallis test
# did i set this up right?
kruskal.test(Environmental$Temp~Environmental$Temperature)# **diff
##
## Kruskal-Wallis rank sum test
##
## data: Environmental$Temp by Environmental$Temperature
## Kruskal-Wallis chi-squared = 203.29, df = 1, p-value < 2.2e-16
kruskal.test(Environmental$pHspec~Environmental$pH)# **diff
##
## Kruskal-Wallis rank sum test
##
## data: Environmental$pHspec by Environmental$pH
## Kruskal-Wallis chi-squared = 470.87, df = 1, p-value < 2.2e-16
kruskal.test(Environmental$pHspec~Environmental$Temperature)# Not dif
##
## Kruskal-Wallis rank sum test
##
## data: Environmental$pHspec by Environmental$Temperature
## Kruskal-Wallis chi-squared = 0.080493, df = 1, p-value = 0.7766
kruskal.test(Environmental$Temp~Environmental$pH)# not dif
##
## Kruskal-Wallis rank sum test
##
## data: Environmental$Temp by Environmental$pH
## Kruskal-Wallis chi-squared = 0.18832, df = 1, p-value = 0.6643
tempsummary<-results %>%
select("Day", "Temperature", "Temp") %>%
drop_na(Temp) %>%
group_by(Day, Temperature) %>%
mutate(meanT = mean(Temp)) %>%
group_by(Temperature) %>%
mutate(sd = sd(Temp)) %>%
mutate(se = sd/sqrt(12))
tempplot<-ggplot(data=tempsummary, aes(Day, meanT, color = Temperature)) +
geom_point(size = 2.5, show.legend=FALSE) +
geom_errorbar(aes(ymin = meanT-se, ymax = meanT+se)) +
geom_line(size=1.2) +
geom_vline(xintercept=0) +
scale_x_continuous(name="Time (Day)", breaks = seq(0,130, by = 10)) +
scale_y_continuous(name="Temperature (°C)", breaks = seq(20,30, by = 1)) +
scale_shape_discrete(name = NULL,
labels = c("Ambient", "High"))+
scale_color_manual(name = NULL,
values = c("gray40", "firebrick3"),
labels = c("Ambient", "High")) +
ggtitle("")+
theme_classic() +
theme(plot.margin=unit(c(0.3,0.6,0.3,0.3), "cm"))+
theme(legend.margin = margin(0),
legend.position = c(.98,.8),
legend.justification = c("right", "top"),
legend.background = element_blank(),
axis.title = element_text(size = 16, face="bold"),
axis.text = element_text(size = 14),
plot.title = element_text(size=16, face="bold"),
legend.text = element_text(size=14, face="bold"));tempplot
#ggsave(filename="Figures/20200611/TempFig1.png", plot=tempplot, dpi=300, width=7, height=5, units="in")
pHsummary<-results %>%
select("Day","pH", "pHspec") %>%
drop_na(pHspec) %>%
group_by(Day, pH) %>%
mutate(meanpH = mean(pHspec)) %>%
group_by(pH, meanpH) %>%
mutate(sd = sd(pHspec)) %>%
mutate(se = sd/sqrt(12))
pHplot<-ggplot(data=pHsummary, aes(Day, meanpH, color = pH)) +
geom_point(size = 2.5, show.legend=FALSE) +
geom_errorbar(aes(ymin = meanpH-sd, ymax = meanpH+sd)) +
geom_line(size=1.2) +
geom_vline(xintercept=0) +
scale_x_continuous(name="Time (Day)", breaks = seq(0,130, by = 10)) +
scale_y_continuous(name="pH (Total Scale)", breaks = seq(7.5,8.1, by = .1)) +
scale_shape_discrete(name = NULL,
labels = c("Low pH", "High pH"))+
scale_color_manual(name = NULL,
values = c("skyblue3", "gray40"),
labels = c("Acidified", "Ambient"),
guide = guide_legend(reverse = TRUE)) +
ggtitle("")+
theme_classic() +
theme(plot.margin=unit(c(0.3,0.6,0.3,0.3), "cm"))+
theme(legend.margin = margin(0),
legend.position = c(.98,.8),
legend.justification = c("right", "top"),
legend.background = element_blank(),
axis.title = element_text(size = 16, face="bold"),
axis.text = element_text(size = 14),
plot.title = element_text(size=16, face="bold"),
legend.text = element_text(size=14, face="bold"));pHplot
#ggsave(filename="Figures/20200611/pHFig1.png", plot=pHplot, dpi=300, width=7, height=5, units="in")
Assemble data for diameters of initial collection of urchins from hatchery and calculate mean
##Initial Collection = Day -24
InitialDiameter <-
#seperate out the initial sizes of urchins on day -24 after initial collection.
results %>%
select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
filter(TankID != "8t" ) %>%
#remove urchin number 8 which died halfway through
drop_na(Diam1, Diam2) %>%
gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
#gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
filter(Day == "-24") %>%
#sizes on day -24
select("TankID", "Diameter", "Temperature", "pH", "Treatment")
InitialMean <- mean(InitialDiameter$Diameter)
InitialSE <- sd(InitialDiameter$Diameter)/sqrt(23)
InitialMean
InitialSE
Assemble Data for diameters of acclimated urchins before conditioning and calculate mean
##After acclimating before ramp up = Day-13
AcclimatedDiameter <-
#seperate out the initial sizes of urchins on day -24, when desired conditions were met.
results %>%
select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
filter(TankID != "8t" ) %>%
#remove urchin number 8 which died halfway through
drop_na(Diam1, Diam2) %>%
gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
#gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
filter(Day == "-13") %>%
#sizes on day -24
select("TankID", "Diameter", "Temperature", "pH", "Treatment")
AcclimatedMean <- mean(AcclimatedDiameter$Diameter)
AcclimatedSE <- sd(AcclimatedDiameter$Diameter)/sqrt(23)
AcclimatedMean
AcclimatedSE
Assemble Data for diameters on day 1 of experiment when desired conditions were reached
##After ramp up and desired conditions are reached to begin experiment = Day 1
Day1Diameter <-
#seperate out the initial sizes of urchins on day 1, when desired conditions were met.
results %>%
select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
filter(TankID != "8t" ) %>%
#remove urchin number 8 which died halfway through
drop_na(Diam1, Diam2) %>%
gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
#gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
filter(Day == "1") %>%
#sizes on day 1
select("TankID", "Diameter", "Temperature", "pH", "Treatment")
Day1Mean <- mean(Day1Diameter$Diameter)
Day1SE <- sd(Day1Diameter$Diameter)/sqrt(23)
Day1Mean
Day1SE
Assemble Data for diameters on day 126 of experiment after full growth
##After ramp up and desired conditions are reached to begin experiment = Day 1
Day126Diameter <-
#seperate out the initial sizes of urchins on day 1, when desired conditions were met.
results %>%
select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
filter(TankID != "8t" ) %>%
#remove urchin number 8 which died halfway through
drop_na(Diam1, Diam2) %>%
gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
#gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
filter(Day == "126") %>%
#sizes on day 1
select("TankID", "Diameter", "Temperature", "pH", "Treatment")
Day126Mean <- mean(Day126Diameter$Diameter)
Day126SE <- sd(Day126Diameter$Diameter)/sqrt(23)
Day126Mean
Day126SE
| Day | Mean Urchin Diameter (mm) | se |
|---|---|---|
| -24 | 7.53 | 0.15 |
| -13 | 10.38 | 0.23 |
| 1 | 16.03 | 0.36 |
| 126 | 70.52 | 1.43 |
Day1Mod <- lmer(Diameter ~ Temperature*pH + (1|TankID), data=Day1Diameter)
anova(Day1Mod, type = 2)
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Temperature 0.04177 0.04177 1 19 0.1922 0.66606
## pH 0.38094 0.38094 1 19 1.7525 0.20128
## Temperature:pH 0.89969 0.89969 1 19 4.1389 0.05611 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of diameter model (Day 1)
#check assumptions
# 1. Normality of residuals
qqPlot(residuals(Day1Mod)) #normal
shapiro.test(residuals(Day1Mod)) #passes
hist(residuals(Day1Mod)) #normal
#check assumptions
# 2. Equal variances
leveneTest(residuals(Day1Mod)~InitialDiameter$Treatment) #Passes
plot(fitted(Day1Mod),resid(Day1Mod,type="pearson"),col="blue")
Assemble data and calculate means using day 1 as initial diameter and calculate means
### GROWTH MODEL: using day 1 as the first official day of experiment when desired environmental conditions were met (24 days after getting urchins. In that 24 days, they were acclimated and then conditions ramped up)
Initial <-
#seperate out the initial sizes of urchins on day 1, when desired conditions were met to standardize growth to this measurement.
results %>%
select("Day", "Treatment", "Temperature", "pH", "Diam1", "Diam2", "TankID") %>%
filter(TankID != "8t" ) %>%
#remove urchin number 8 which died halfway through
drop_na(Diam1, Diam2) %>%
gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
#gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously received)
filter(Day == "1") %>%
#sizes on day 1
select("Temperature", "pH", "Diameter", "TankID")%>%
group_by(TankID)%>%
mutate(Diameter=mean(Diameter))%>%
rename(InitialDiameter=Diameter)%>%
unique()
Final <-
#seperate out the final sizes of urchins on day 126
results %>%
select("Day", "Treatment", "Temperature", "pH", "Diam1", "Diam2", "TankID") %>%
drop_na(Diam1, Diam2) %>%
gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
#gather Diam1 and Diam2 into single column
filter(Day == "126") %>%
#filter for sizes on day 126
select("Temperature", "pH", "Diameter","TankID")%>%
group_by(TankID)%>%
mutate(Diameter=mean(Diameter))%>%
rename(FinalDiameter=Diameter)%>%
unique()
resultsGrow <-
#create column that is growth
full_join(Initial, Final) %>%
#binds initial (Diameter) and final (Diameter1) sizes to calculate growth
mutate(Growth = ((FinalDiameter - InitialDiameter)/InitialDiameter*100)) %>%
#add column for growth calculation
select("Temperature", "pH", "TankID", "Growth")
meanGrowth <-
resultsGrow %>%
dplyr::group_by(Temperature, pH) %>%
dplyr::summarise(meanGrowth = mean(Growth), s.e. = se(Growth), n=length(Growth))
#Figure out means of each treatment group
| Temp | pH | Mean Growth (%) | se |
|---|---|---|---|
| Amb | OA | 326.62 | 0.15 |
| Amb | Amb | 323.25 | 0.23 |
| Heat | OA | 352.02 | 0.36 |
| Heat | Amb | 376.25 | 1.43 |
FINAL growth ANOVA results
#LMM for growth:
modGrow <-
resultsGrow %>%
aov(Growth~ Temperature * pH, data=.)
summary(modGrow)
## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 1 8375 8375 4.038 0.0589 .
## pH 1 545 545 0.263 0.6141
## Temperature:pH 1 1082 1082 0.522 0.4789
## Residuals 19 39403 2074
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions for FINAL growth model
## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modGrow))
## [1] 15 18
shapiro.test(residuals(modGrow)) #passes
##
## Shapiro-Wilk normality test
##
## data: residuals(modGrow)
## W = 0.96259, p-value = 0.5175
hist(residuals(modGrow)) #looks normal
## ASSSUMPTIONS
# 2. Equal variances
leveneTest(residuals(modGrow)~resultsGrow$Temperature) #Passes
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.2731 0.6068
## 21
leveneTest(residuals(modGrow)~resultsGrow$pH) #doesn't pass
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 3.9714 0.05944 .
## 21
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(fitted(modGrow),resid(modGrow,type="pearson"),col="blue")
Power Analysis for FINAL growth model: Calculate Effect Size
#calculate effect size for difference in growth by temperature.
library("effsize")
library("rstatix")
#calculate cohens D effect size
#cohens_d(Growth ~ Temperature, var.equal = FALSE, ref.group="ambient", data=resultsGrow)
#somethings acting weird here, but now not using this analysis
#From this, we see that we have a large effect size between temperature treatments ().
Power Analysis for FINAL growth model: Confirm Sample size
#Confirm that we have the sapmle size necessary to detect a large effect. Power analysis for test
#calculate test power first given known variables
pwr.anova.test(k=4, n=6, f=0.8, sig.level = 0.05) #balanced one way anova
##
## Balanced one-way analysis of variance power calculation
##
## k = 4
## n = 6
## f = 0.8
## sig.level = 0.05
## power = 0.8592125
##
## NOTE: n is number in each group
#k=Number of groups - 4 treatments
#n=Number of observations (per group) - 6
#f=Effect size - 0.8 based on above calculation
#sig.level=Significance level (Type I error probability)
#solve for power=Power of test (1 minus Type II error probability)
#this last test identified power as 0.86, now look for the sample size we would need to detect a small (0.2), medium (0.5) or large (0.8) effect size.
pwr.anova.test(k=4, f=0.2, sig.level = 0.05, power=0.86) #n needed for small effect, would need n=79 urchins
##
## Balanced one-way analysis of variance power calculation
##
## k = 4
## n = 79.90233
## f = 0.2
## sig.level = 0.05
## power = 0.86
##
## NOTE: n is number in each group
pwr.anova.test(k=4, f=0.5, sig.level = 0.05, power=0.86) #n needed for small effect, would need n=13 urchins
##
## Balanced one-way analysis of variance power calculation
##
## k = 4
## n = 13.64663
## f = 0.5
## sig.level = 0.05
## power = 0.86
##
## NOTE: n is number in each group
pwr.anova.test(k=4, f=0.8, sig.level = 0.05, power=0.86) #n needed for small effect, would need n=6 urchins
##
## Balanced one-way analysis of variance power calculation
##
## k = 4
## n = 6.010194
## f = 0.8
## sig.level = 0.05
## power = 0.86
##
## NOTE: n is number in each group
#This power analysis for growth shows us that in order to detect a small effect, we would need 79 urchins. To detect a medium effect, we would need 13 urchins. To detect a large effect, we would need 6 urchins. Given that we have n=6 urchins per group, we would have the statistical power to detect a large effect that was detected. This validates our test power.
growthsummary<-resultsGrow %>%
#filter(Day==126) %>%
select("pH", "Temperature","Growth") %>%
drop_na(Growth) %>%
mutate(meanPct = Growth) %>%
group_by(pH, Temperature) %>%
summarise(mean = mean(meanPct),
N = length(meanPct),
se = std.error(meanPct))
growthplot<-ggplot(data=growthsummary, aes(x=Temperature, y=mean, color = pH, group=interaction(Temperature,pH))) +
geom_point(size = 3, position=position_dodge(0.1)) +
geom_line(aes(group=pH), position=position_dodge(0.1)) +
xlab("Temperature Treatment") +
ylab("Total Growth (%)")+
geom_errorbar(aes(ymin = mean-se, ymax = mean+se, width = 0.1), position=position_dodge(0.1)) +
theme_bw() +
scale_x_discrete(limits=c("ambient", "heated"),
label=c("Ambient", "High"))+
scale_colour_manual(name = "pH Treatment",
values = c("skyblue3", "gray40"),
labels = c("Acidified", "Ambient"),
guide = guide_legend(reverse = TRUE)) +
#geom_text(x=1.5, y=125, label="p(pH)=0.466", size=6, color="darkgray") +
#geom_text(x=1.5, y=100, label="p(Temperature)=0.006", size=6, color="black") +
#geom_text(x=1.5, y=75, label="p(Interaction)=0.303", size=6, color="darkgray") +
ylim(0,450)+
theme_classic()+
theme(legend.position = "right",
legend.background = element_blank(),
axis.title = element_text(size = 16, face="bold"),
axis.text = element_text(size = 14),
plot.title = element_text(size=16, face="bold"),
legend.text = element_text(size=14),
legend.title = element_text(size=16, face="bold"));growthplot
#ggsave(filename="Figures/20200611/GrowthFig.png", plot=growthplot, dpi=300, width=8, height=6, units="in")
resultsGrowTime <- #create table of growth in treatments
results %>%
select("Day", "TankID", "Treatment","Temperature","pH","Growth") %>%
drop_na(Growth)
#doing square transformation here to account for smaller values at the start of the experiment
modelGrowTime <- lmer(Growth^2~Temperature*pH*Day + (1|TankID), data=resultsGrowTime)
summary(modelGrowTime) #generate results
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Growth^2 ~ Temperature * pH * Day + (1 | TankID)
## Data: resultsGrowTime
##
## REML criterion at convergence: 1284.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2189 -0.6094 -0.1552 0.5379 4.4363
##
## Random effects:
## Groups Name Variance Std.Dev.
## TankID (Intercept) 0.9657 0.9827
## Residual 1.3313 1.1538
## Number of obs: 382, groups: TankID, 24
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -1.215915 0.456538 28.351183 -2.663
## Temperatureheated -0.409572 0.645642 28.351183 -0.634
## pHambient -0.075016 0.645642 28.351183 -0.116
## Day 0.093002 0.003091 354.008336 30.093
## Temperatureheated:pHambient -0.231471 0.913492 28.401729 -0.253
## Temperatureheated:Day 0.015997 0.004371 354.008336 3.660
## pHambient:Day -0.001069 0.004371 354.008336 -0.244
## Temperatureheated:pHambient:Day 0.010159 0.006230 354.116662 1.631
## Pr(>|t|)
## (Intercept) 0.012615 *
## Temperatureheated 0.530929
## pHambient 0.908322
## Day < 2e-16 ***
## Temperatureheated:pHambient 0.801791
## Temperatureheated:Day 0.000291 ***
## pHambient:Day 0.806997
## Temperatureheated:pHambient:Day 0.103844
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tmprtr pHmbnt Day Tmpr:H Tmpr:D pHmb:D
## Tempertrhtd -0.707
## pHambient -0.707 0.500
## Day -0.402 0.284 0.284
## Tmprtrhtd:H 0.500 -0.707 -0.707 -0.201
## Tmprtrhtd:D 0.284 -0.402 -0.201 -0.707 0.284
## pHambint:Dy 0.284 -0.201 -0.402 -0.707 0.284 0.500
## Tmprtrh:H:D -0.199 0.282 0.282 0.496 -0.402 -0.702 -0.702
anova(modelGrowTime, type=2)
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Temperature 1.8 1.8 1 28.40 1.3206 0.2601
## pH 0.2 0.2 1 28.40 0.1737 0.6800
## Day 5804.6 5804.6 1 354.11 4360.0148 < 2.2e-16 ***
## Temperature:pH 0.1 0.1 1 28.40 0.0642 0.8018
## Temperature:Day 60.6 60.6 1 354.11 45.4855 6.275e-11 ***
## pH:Day 2.1 2.1 1 354.11 1.6002 0.2067
## Temperature:pH:Day 3.5 3.5 1 354.12 2.6592 0.1038
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of model.
hist(residuals(modelGrowTime))
qqPlot(residuals(modelGrowTime))
## [1] 377 374
resultsGrowTime$code<-paste(resultsGrowTime$Treatment, resultsGrowTime$Day)
leveneTest(residuals(modelGrowTime)~resultsGrowTime$code)
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 63 2.1709 6.864e-06 ***
## 318
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Power Analysis for growth OVER TIME model: Calculate Effect Size
#calculate effect size for difference in growth by temperature.
library("effsize")
library("rstatix")
#calculate cohens D effect size
cohens_d(Growth ~ Temperature, var.equal = FALSE, ref.group="ambient", data=resultsGrowTime)
## # A tibble: 1 × 7
## .y. group1 group2 effsize n1 n2 magnitude
## * <chr> <chr> <chr> <dbl> <int> <int> <ord>
## 1 Growth ambient heated -0.100 192 190 negligible
#From this, we see that we have a negligible effect size between temperature treatments (-0.100).
Power Analysis for OVER TIME growth model: Confirm Sample size
#Confirm that we have the sample size necessary to detect a large effect. Power analysis for test
#calculate test power first given known variables
pwr.anova.test(k=4, n=96, f=0.1, sig.level = 0.05) #balanced one way anova
##
## Balanced one-way analysis of variance power calculation
##
## k = 4
## n = 96
## f = 0.1
## sig.level = 0.05
## power = 0.3418767
##
## NOTE: n is number in each group
#k=Number of groups - 4 treatments
#n=Number of observations (per group) - 96
resultsGrowTime %>%
count(Treatment)
## Treatment n
## 1 AMB 96
## 2 OA 96
## 3 T 94
## 4 T/OA 96
#f=Effect size - 0.1 based on above calculation
#sig.level=Significance level (Type I error probability)
#solve for power=Power of test (1 minus Type II error probability)
#this last test identified power as 0.34, now look for the sample size we would need to detect a small (0.2), medium (0.5) or large (0.8) effect size.
pwr.anova.test(k=4, f=0.2, sig.level = 0.05, power=0.34) #n needed for small effect, would need n=25 urchins
##
## Balanced one-way analysis of variance power calculation
##
## k = 4
## n = 24.60756
## f = 0.2
## sig.level = 0.05
## power = 0.34
##
## NOTE: n is number in each group
pwr.anova.test(k=4, f=0.5, sig.level = 0.05, power=0.34) #n needed for small effect, would need n=5 urchins
##
## Balanced one-way analysis of variance power calculation
##
## k = 4
## n = 4.839745
## f = 0.5
## sig.level = 0.05
## power = 0.34
##
## NOTE: n is number in each group
pwr.anova.test(k=4, f=0.8, sig.level = 0.05, power=0.34) #n needed for small effect, would need n=3 urchins
##
## Balanced one-way analysis of variance power calculation
##
## k = 4
## n = 2.614108
## f = 0.8
## sig.level = 0.05
## power = 0.34
##
## NOTE: n is number in each group
#This power analysis for growth over time shows us that in order to detect a small effect, we would need 25 urchins. To detect a medium effect, we would need 5 urchins. To detect a large effect, we would need 3 urchins. Given that we have n=6 urchins per group, we would have the statistical power to detect a medium effect that was detected. This validates our test power.
Figure 4: Growth over time across all treatments
#### Growth Across Treatments
results %>%
select("Day", "Treatment","Temperature","pH","Growth") %>%
drop_na(Growth) %>%
group_by(Day, Treatment) %>%
summarise(se = (std.error(Growth))*100, meanPct = (mean(Growth)*100)) %>%
ggplot(aes(Day, meanPct, shape = Treatment, color = Treatment)) +
geom_point(size = 2.5) +
geom_line() +
scale_x_continuous(name="Time (Day)", breaks = seq(0,130, by = 10)) +
scale_y_continuous(name="Growth (%)", breaks = seq(0,400, by = 25)) +
geom_errorbar(aes(ymin = meanPct-se, ymax = meanPct+se, width = 1.5)) +
theme_classic() +
scale_shape_discrete(name = NULL,
labels = c("Amb T, Amb pH", "Amb T, Low pH", "High T, Amb pH", "High T, Low pH")) +
scale_colour_manual(name = NULL,
values = c("gray40", "skyblue3", "firebrick3","mediumpurple3"),
labels = c("Amb T, Amb pH", "Amb T, Low pH", "High T, Amb pH", "High T, Low pH")) +
theme(legend.margin = margin(0),
legend.position = c(.1, .98),
legend.justification = c("left", "top"),
legend.background = element_blank())
#ggsave(filename="Figures/20200611/GrowthTreatments.png", dpi=300, width=8, height=6, units="in")
Figure 5: Growth over time plot by temperature
GrowthTemp <-
results %>%
select("Day","Temperature","Growth") %>%
drop_na(Growth) %>%
group_by(Day, Temperature) %>%
summarise(se = (std.error(Growth))*100, meanPct = (mean(Growth)*100)) %>%
ggplot(aes(Day, meanPct, shape = Temperature, color = Temperature)) +
geom_point(size = 2.5) +
geom_line() +
scale_x_continuous(name="Time (Day)", breaks = seq(0,130, by = 10)) +
scale_y_continuous(name="Growth (%)", breaks = seq(0,400, by = 25)) +
geom_errorbar(aes(ymin = meanPct-se, ymax = meanPct+se, width = 1.5)) +
theme_classic() +
scale_colour_manual(name = NULL,
values = c("gray40", "firebrick3"),
labels = c("Ambient Temperture", "High Temperature")) +
scale_shape_discrete(name = NULL,
labels = c("Ambient Temperture", "High Temperature")) +
theme(legend.margin = margin(0),
legend.position = c(.03, .93),
legend.justification = c("left", "top")) +
annotate("text", x = 0, y = 375, label = "(a)") +
annotate("text", x = 130, y = 359, label = "*", size = 7); GrowthTemp
## `summarise()` has grouped output by 'Day'. You can override using the `.groups` argument.
#ggsave(filename="Figures/20200611/GrowthTemp.png", plot=GrowthTemp, dpi=300, width=8, height=6, units="in")
Figure 6. Growth over time plot by pH
GrowthpH <-
results %>%
select("Day","pH","Growth") %>%
drop_na(Growth) %>%
group_by(Day, pH) %>%
summarise(se = (std.error(Growth))*100, meanPct = (mean(Growth)*100)) %>%
ggplot(aes(Day, meanPct, shape = pH, color = pH)) +
geom_point(size = 2.5) +
geom_line() +
scale_x_continuous(name="Time (Day)", breaks = seq(0,130, by = 10)) +
scale_y_continuous(name="Growth (%)", breaks = seq(0,400, by = 25)) +
geom_errorbar(aes(ymin = meanPct-se, ymax = meanPct+se, width = 1.5)) +
theme_classic() +
scale_colour_manual(name = NULL,
values = c("skyblue3", "gray40"),
labels = c("Low pH", "Ambient pH")) +
scale_shape_discrete(name = NULL,
labels = c("Low pH", "Ambient pH")) +
theme(legend.margin = margin(0),
legend.position = c(.03, .93),
legend.justification = c("left", "top")) +
annotate("text", x = 0, y = 370, label = "(b)"); GrowthpH
## `summarise()` has grouped output by 'Day'. You can override using the `.groups` argument.
ggsave(filename="Figures/20200611/GrowthpH.png", plot=GrowthpH, dpi=300, width=8, height=6, units="in")
Assemble data for chosen spine tips (distal end)
### Model 1: calcification at the tips of spines
resultsTip <-
#create data using only the SEM images of the tips of spines.
results %>%
select("Treatment", "Temperature", "pH", "PartOfSpine", "Chosen", "RatioSEM", "TankID") %>%
#selects for desired variables. PartOfSpine is tip or base, Chosen indicates which side of the image (arbitrary) is selected for (aka has no dust)
drop_na(RatioSEM) %>%
filter(Chosen == "yes", PartOfSpine == "Tip")
#filters only for tips and chosen side. No dusty images to be analyzed.
Results for spine tips with linear mixed model
#LMM for calcification at the tips of spines
modTip <-
resultsTip %>%
lmer(RatioSEM ~ Temperature * pH + (1|TankID), data = .)
#somethign in this model is weird - treating each measure as individual and not as part of the random effect?
summary(modTip)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RatioSEM ~ Temperature * pH + (1 | TankID)
## Data: .
##
## REML criterion at convergence: 64.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.76413 -0.74507 -0.02052 0.63045 2.07630
##
## Random effects:
## Groups Name Variance Std.Dev.
## TankID (Intercept) 0.0000 0.0000
## Residual 0.1356 0.3682
## Number of obs: 67, groups: TankID, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.69118 0.08930 63.00000 18.938 <2e-16 ***
## Temperatureheated -0.06662 0.12453 63.00000 -0.535 0.5945
## pHambient -0.07462 0.12453 63.00000 -0.599 0.5512
## Temperatureheated:pHambient 0.30557 0.18089 63.00000 1.689 0.0961 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.717
## pHambient -0.717 0.514
## Tmprtrhtd:H 0.494 -0.688 -0.688
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
anova(modTip, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type II Analysis of Variance Table with Satterthwaite's method.
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Temperature 0.10158 0.10158 1 63 0.7492 0.39000
## pH 0.08185 0.08185 1 63 0.6038 0.44006
## Temperature:pH 0.38685 0.38685 1 63 2.8534 0.09612 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions for calcification at tip model.
## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modTip)) #Normal
## [1] 38 8
shapiro.test(residuals(modTip)) #Pass
##
## Shapiro-Wilk normality test
##
## data: residuals(modTip)
## W = 0.97737, p-value = 0.2605
hist(residuals(modTip))
# 2. Equal variances
leveneTest(residuals(modTip)~resultsTip$Treatment) #Pass
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 0.4545 0.715
## 63
plot(fitted(modTip),resid(modTip,type="pearson"),col="blue")
Power tests: Effect Size
library("effsize")
library("rstatix")
#calculate cohens D effect size
cohens_d(RatioSEM ~ pH, var.equal = FALSE, ref.group="ambient", data=resultsTip)
## # A tibble: 1 × 7
## .y. group1 group2 effsize n1 n2 magnitude
## * <chr> <chr> <chr> <dbl> <int> <int> <ord>
## 1 RatioSEM ambient acidified 0.172 32 35 negligible
#From this, we see that we have a small effect size between pH treatments (0.17).
Power test: Confirm sample size
#calculate test power first given known variables
pwr.anova.test(k=4, n=18, f=0.17, sig.level = 0.05) #balanced one way anova
##
## Balanced one-way analysis of variance power calculation
##
## k = 4
## n = 18
## f = 0.17
## sig.level = 0.05
## power = 0.1894901
##
## NOTE: n is number in each group
#k=Number of groups - 4 treatments
#n=Number of observations (per group) - 18 tips per treatment group
#f=Effect size - 0.17 based on above calculation
#sig.level=Significance level (Type I error probability)
#solve for power=Power of test (1 minus Type II error probability)
#this last test identified power as 0.19, now look for the sample size we would need to detect a small (0.2), medium (0.5) or large (0.8) effect size.
pwr.anova.test(k=4, f=0.2, sig.level = 0.05, power=0.19) #n needed for small effect, would need n=13 spines
##
## Balanced one-way analysis of variance power calculation
##
## k = 4
## n = 13.32761
## f = 0.2
## sig.level = 0.05
## power = 0.19
##
## NOTE: n is number in each group
pwr.anova.test(k=4, f=0.5, sig.level = 0.05, power=0.19) #n needed for medium effect, would need n=3 spines
##
## Balanced one-way analysis of variance power calculation
##
## k = 4
## n = 3.061966
## f = 0.5
## sig.level = 0.05
## power = 0.19
##
## NOTE: n is number in each group
#pwr.anova.test(k=4, f=0.8, sig.level = 0.05, power=0.19) #n needed for large effect, would need n=2 spines - this one doesn't work, maybe showing you don't need any urchins?
Assemble data for chosen spine bases (proximal end)
### Model 2: calcification in the base of the spines
resultsBase <-
#create data using only the SEM images of the base of spines.
results %>%
select("Treatment", "Temperature", "pH", "PartOfSpine", "Chosen", "RatioSEM", "TankID") %>%
#selects for desired variables. PartOfSpine is tip or base, Chosen indicates which side of the image (arbitrary) is selected for (aka has no dust)
drop_na(RatioSEM) %>%
filter(Chosen == "yes", PartOfSpine == "Base")
#filters only for tips and chosen side. No dusty images to be analyzed.
Results for spine bases with linear mixed model
#LMM for calcification at the tips of spines
modBase <-
resultsBase %>%
lmer(RatioSEM~ Temperature * pH + (1|TankID), data=.)
summary(modBase)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RatioSEM ~ Temperature * pH + (1 | TankID)
## Data: .
##
## REML criterion at convergence: 98.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5535 -0.6315 -0.1529 0.6918 1.9937
##
## Random effects:
## Groups Name Variance Std.Dev.
## TankID (Intercept) 0.1303 0.3610
## Residual 0.1597 0.3996
## Number of obs: 68, groups: TankID, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.1440 0.1749 18.8077 12.257 2.06e-10 ***
## Temperatureheated 0.3566 0.2487 19.1707 1.434 0.1677
## pHambient 0.6524 0.2474 18.8077 2.637 0.0163 *
## Temperatureheated:pHambient -0.2224 0.3594 18.9804 -0.619 0.5434
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.703
## pHambient -0.707 0.497
## Tmprtrhtd:H 0.487 -0.692 -0.688
anova(modBase, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type III Analysis of Variance Table with Satterthwaite's method.
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Temperature 0.30743 0.30743 1 18.992 1.9251 0.181362
## pH 1.48873 1.48873 1 18.960 9.3223 0.006553 **
## Temperature:pH 0.06114 0.06114 1 18.980 0.3829 0.543424
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions for calcification at base model.
## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modBase)) #Normal
## [1] 33 37
shapiro.test(residuals(modBase)) # Pass
##
## Shapiro-Wilk normality test
##
## data: residuals(modBase)
## W = 0.97304, p-value = 0.1472
hist(residuals(modBase))
# 2. Equal variances
leveneTest(residuals(modBase)~resultsBase$Treatment) #Pass
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 1.068 0.369
## 64
plot(fitted(modBase),resid(modBase,type="pearson"),col="blue")
Check power of the test
library("effsize")
library("rstatix")
#calculate cohens D effect size
cohens_d(RatioSEM ~ pH, var.equal = FALSE, ref.group="ambient", data=resultsBase)
## # A tibble: 1 × 7
## .y. group1 group2 effsize n1 n2 magnitude
## * <chr> <chr> <chr> <dbl> <int> <int> <ord>
## 1 RatioSEM ambient acidified 1.03 33 35 large
From this, we see that we have a large effect size between temperature treatments (1.03).
Confirm that we have the sample size necessary to detect a large effect. Power analysis for test
#calculate test power first given known variables
pwr.anova.test(k=4, n=18, f=1.03, sig.level = 0.05) #balanced one way anova
##
## Balanced one-way analysis of variance power calculation
##
## k = 4
## n = 18
## f = 1.03
## sig.level = 0.05
## power = 1
##
## NOTE: n is number in each group
#k=Number of groups - 4 treatments
#n=Number of observations (per group) - 18 bases per treatment group
#f=Effect size - 0.20 based on above calculation
#sig.level=Significance level (Type I error probability)
#solve for power=Power of test (1 minus Type II error probability)
#this last test identified power as 1, now look for the sample size we would need to detect a small (0.2), medium (0.5) or large (0.8) effect size.
pwr.anova.test(k=4, f=0.2, sig.level = 0.05, power=1) #n needed for small effect, would need n= spines
##
## Balanced one-way analysis of variance power calculation
##
## k = 4
## n = 1e+09
## f = 0.2
## sig.level = 0.05
## power = 1
##
## NOTE: n is number in each group
pwr.anova.test(k=4, f=0.5, sig.level = 0.05, power=1) #n needed for small effect, would need n= spines
##
## Balanced one-way analysis of variance power calculation
##
## k = 4
## n = 1e+09
## f = 0.5
## sig.level = 0.05
## power = 1
##
## NOTE: n is number in each group
pwr.anova.test(k=4, f=0.8, sig.level = 0.05, power=1) #n needed for small effect, would need n= spines
##
## Balanced one-way analysis of variance power calculation
##
## k = 4
## n = 1e+09
## f = 0.8
## sig.level = 0.05
## power = 1
##
## NOTE: n is number in each group
Make Interaction plots for calcification at the base and the tip
#### Tip and base calcification ratios
#base
basesummary<-results %>%
select("Temperature","pH","RatioSEM", "PartOfSpine") %>%
filter(PartOfSpine=="Base")%>%
drop_na(RatioSEM) %>%
group_by(pH, Temperature, PartOfSpine) %>%
summarize(mean = mean(RatioSEM), se = std.error(RatioSEM), N = length(RatioSEM))
baseplot<- ggplot(data=basesummary, aes(x=Temperature, y=mean, color = pH, group=interaction(Temperature,pH))) +
geom_point(size = 3, position=position_dodge(0.1)) +
geom_line(aes(group=pH), position=position_dodge(0.1)) +
xlab("Temperature Treatment") +
ylab("Proximal Calcification Ratio")+
ylim(1,3.5)+
geom_errorbar(aes(ymin = mean-se, ymax = mean+se, width = 0.1), position=position_dodge(0.1)) +
theme_bw() +
scale_x_discrete(limits=c("ambient", "heated"),
label=c("Ambient", "High"))+
scale_colour_manual(name = "pH Treatment",
values = c("skyblue3", "gray40"),
labels = c("Acidified", "Ambient"),
guide = guide_legend(reverse = TRUE)) +
geom_text(x=1.5, y=1.6, label="p(pH)=0.006", size=6, color="black") +
geom_text(x=1.5, y=1.4, label="p(Temperature)=0.18", size=6, color="darkgray") +
geom_text(x=1.5, y=1.2, label="p(Interaction)=0.54", size=6, color="darkgray") +
theme_classic()+
theme(legend.position = "right",
legend.background = element_blank(),
axis.title = element_text(size = 16, face="bold"),
axis.text = element_text(size = 14),
plot.title = element_text(size=16, face="bold"),
legend.text = element_text(size=14),
legend.title = element_text(size=16, face="bold"));baseplot
ggsave(filename="Figures/20200611/BaseCalcificationFig.png", plot=baseplot, dpi=300, width=8, height=6, units="in")
#tip
tipsummary<-results %>%
select("Temperature","pH","RatioSEM", "PartOfSpine") %>%
filter(PartOfSpine=="Tip")%>%
drop_na(RatioSEM) %>%
group_by(pH, Temperature, PartOfSpine) %>%
summarize(mean = mean(RatioSEM), se = std.error(RatioSEM), N = length(RatioSEM))
tipplot<-ggplot(data=tipsummary, aes(x=Temperature, y=mean, color = pH, group=interaction(Temperature,pH))) +
geom_point(size = 3, position=position_dodge(0.1)) +
geom_line(aes(group=pH), position=position_dodge(0.1)) +
xlab("Temperature Treatment") +
ylab("Distal Calcification Ratio")+
ylim(1,3.5)+
geom_errorbar(aes(ymin = mean-se, ymax = mean+se, width = 0.1), position=position_dodge(0.1)) +
theme_bw() +
scale_x_discrete(limits=c("ambient", "heated"),
label=c("Ambient", "High"))+
scale_colour_manual(name = "pH Treatment",
values = c("skyblue3", "gray40"),
labels = c("Acidified", "Ambient"),
guide = guide_legend(reverse = TRUE)) +
geom_text(x=1.5, y=3.0, label="p(pH)=0.44", size=6, color="darkgray") +
geom_text(x=1.5, y=2.8, label="p(Temperature)=0.39", size=6, color="darkgray") +
geom_text(x=1.5, y=2.6, label="p(Interaction)=0.10", size=6, color="darkgray") +
theme_classic()+
theme(legend.position = "none", #removed the legend for the tip, so when we align them horizontally we only have one legend. But need to adjust the size so it is the same as base
legend.background = element_blank(),
axis.title = element_text(size = 16, face="bold"),
axis.text = element_text(size = 14),
plot.title = element_text(size=16, face="bold"),
legend.text = element_text(size=14),
legend.title = element_text(size=16, face="bold"));tipplot
ggsave(filename="Figures/20200611/TipCalcificationFig.png", plot=tipplot, dpi=300, width=8, height=6, units="in")
Combine both interaction plots at base and tip to be in horizontal line
Assemble data for dropped spines
## Number of loose spines counted at the bottom of tanks at the end of the experiment. These were either broken or shed. Note: these spines were not collected, just counted, so can not determine whether they were broken or shed.
resultsDropped <-
#create data using only the needed variables.
results %>%
select ("Treatment", "Temperature", "pH", "SpineCount", "TankID") %>%
drop_na(SpineCount) %>%
group_by(pH) %>%
mutate(mean = mean(SpineCount), se = std.error(SpineCount), N = length(SpineCount))
Analyse dropped spines with linear mixed effect model.
##LM of dropped spines - can't make it LMM because error: number of levels of each grouping factor must be < number of observations - only have one count for each TankID.
modDropped <-
resultsDropped %>%
lm(SpineCount~ Temperature * pH, data=.)
summary(modDropped)
##
## Call:
## lm(formula = SpineCount ~ Temperature * pH, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.0000 -6.0833 -0.1667 0.4000 20.0000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.167 3.640 5.814 1.33e-05 ***
## Temperatureheated -1.167 5.148 -0.227 0.82315
## pHambient -21.000 5.148 -4.079 0.00064 ***
## Temperatureheated:pHambient 1.600 7.461 0.214 0.83248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.917 on 19 degrees of freedom
## Multiple R-squared: 0.6088, Adjusted R-squared: 0.547
## F-statistic: 9.855 on 3 and 19 DF, p-value: 0.0003899
anova(modDropped) #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type III Analysis of Variance Table with Satterthwaite's method.
## Analysis of Variance Table
##
## Response: SpineCount
## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 1 1.52 1.52 0.0192 0.8914
## pH 1 2345.78 2345.78 29.4995 3.062e-05 ***
## Temperature:pH 1 3.66 3.66 0.0460 0.8325
## Residuals 19 1510.87 79.52
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions for model for dropped spines
##ASSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modDropped)) #not normal
## [1] 2 4
shapiro.test(residuals(modDropped)) #faaaiiilll
##
## Shapiro-Wilk normality test
##
## data: residuals(modDropped)
## W = 0.86706, p-value = 0.005636
hist(residuals(modDropped))
# 2. Equal variances
leveneTest(residuals(modDropped)~resultsDropped$Treatment) # barely pass
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 2.8487 0.06483 .
## 19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(fitted(modDropped),resid(modDropped,type="pearson"),col="blue")
Analysis for dropped spines does not meet assumptions. Attempt a square root transformation.
##DOES NOT MEET ASSUMPTIONS so ATTEMPT AT Transformation:
resultsDropped$tdata<-(resultsDropped$SpineCount)^1/2
# transformation
modDropped1 <- lm(tdata~ Temperature * pH, data=resultsDropped)
# run lm model of transformed data
###ASSUMPTIONS for Transformation
# 1. Normality of residuals
qqPlot(residuals(modDropped1)) #not normal
## [1] 2 4
shapiro.test(residuals(modDropped1)) #fail
##
## Shapiro-Wilk normality test
##
## data: residuals(modDropped1)
## W = 0.86706, p-value = 0.005636
hist(residuals(modDropped1))
# 2. equal variances
bartlett.test(residuals(modDropped1)~resultsDropped$Treatment) #fail
##
## Bartlett test of homogeneity of variances
##
## data: residuals(modDropped1) by resultsDropped$Treatment
## Bartlett's K-squared = 43.282, df = 3, p-value = 2.144e-09
plot(fitted(modDropped1),resid(modDropped1,type="pearson"),col="blue")
summary(modDropped1)
##
## Call:
## lm(formula = tdata ~ Temperature * pH, data = resultsDropped)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.0000 -3.0417 -0.0833 0.2000 10.0000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.5833 1.8202 5.814 1.33e-05 ***
## Temperatureheated -0.5833 2.5742 -0.227 0.82315
## pHambient -10.5000 2.5742 -4.079 0.00064 ***
## Temperatureheated:pHambient 0.8000 3.7304 0.214 0.83248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.459 on 19 degrees of freedom
## Multiple R-squared: 0.6088, Adjusted R-squared: 0.547
## F-statistic: 9.855 on 3 and 19 DF, p-value: 0.0003899
Anova(modDropped1, type=2) #test for significance of model
## Anova Table (Type II tests)
##
## Response: tdata
## Sum Sq Df F value Pr(>F)
## Temperature 0.23 1 0.0118 0.9146
## pH 586.44 1 29.4995 3.062e-05 ***
## Temperature:pH 0.91 1 0.0460 0.8325
## Residuals 377.72 19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Transfomraiton (of any usual kind) of dropped spines was not successful, so we use a non parametric Kruskal Wallis Test.
### Transformations not work, so use non parametric:
kruskal.test(resultsDropped$SpineCount~resultsDropped$Treatment)# SIG*** (p=0.0005563)
kruskal.test(resultsDropped$SpineCount~resultsDropped$Temperature)#
kruskal.test(resultsDropped$SpineCount~resultsDropped$pH)#
Conduct Dunn Post Hoc Test for dropped spines
## POST HOC Pairwise
dunn <-
resultsDropped %>%
dunnTest(SpineCount ~ Treatment,
method = "holm", kw = TRUE, data = .)
dunnph <- dunn$res
cldList(P.adj ~ Comparison, data = dunnph, threshold = 0.05)
## Group Letter MonoLetter
## 1 AMB a a
## 2 OA b b
## 3 T ac a c
## 4 T/OA bc bc
droppedsummary<-results %>%
select("Temperature","pH","SpineCount") %>%
drop_na(SpineCount) %>%
group_by(pH, Temperature) %>%
summarize(mean = mean(SpineCount), se = std.error(SpineCount), N = length(SpineCount))
dropplot<-ggplot(data=droppedsummary, aes(x=Temperature, y=mean, color = pH, group=interaction(Temperature,pH))) +
geom_point(size = 3, position=position_dodge(0.1)) +
geom_line(aes(group=pH), position=position_dodge(0.1)) +
xlab("Temperature Treatment") +
ylab("Dropped Spines")+
geom_errorbar(aes(ymin = mean-se, ymax = mean+se, width = 0.1), position=position_dodge(0.1)) +
ylim(0,30) +
theme_bw() +
scale_x_discrete(limits=c("ambient", "heated"),
label=c("Ambient", "High"))+
scale_colour_manual(name = "pH Treatment",
values = c("skyblue3", "gray40"),
labels = c("Acidified", "Ambient"),
guide = guide_legend(reverse = TRUE)) +
#geom_text(x=1.5, y=11, label="p(pH)<0.001", size=6, color="black") +
#geom_text(x=1.5, y=9, label="p(Temperature)=0.91", size=6, color="darkgray") +
#geom_text(x=1.5, y=7, label="p(Interaction)=0.83", size=6, color="darkgray") +
theme_classic()+
theme(legend.position = "right",
legend.background = element_blank(),
axis.title = element_text(size = 16, face="bold"),
axis.text = element_text(size = 14),
plot.title = element_text(size=16, face="bold"),
legend.text = element_text(size=14),
legend.title = element_text(size=16, face="bold"));dropplot
#ggsave(filename="Figures/20200611/DroppedSpinesFig.png", plot=dropplot, dpi=300, width=8, height=6, units="in")
#library("effsize")
#library("rstatix")
#calculate cohens D effect size
#cohens_d(SpineCount ~ pH, var.equal = FALSE, ref.group="ambient", data=resultsDropped)
From this, we see that we have a small effect size between temperature treatments (-2.435).
Confirm that we have the sample size necessary to detect a large effect. Power analysis for test
#calculate test power first given known variables
##pwr.anova.test(k=4, n=18, f=2.44, sig.level = 0.05) #balanced one way anova
#k=Number of groups - 4 treatments
#n=Number of observations (per group) - 18 tips per treatment group
#f=Effect size - 0.20 based on above calculation
#sig.level=Significance level (Type I error probability)
#solve for power=Power of test (1 minus Type II error probability)
#this last test identified power as 0.25, now look for the sample size we would need to detect a small (0.2), medium (0.5) or large (0.8) effect size.
##pwr.anova.test(k=4, f=0.2, sig.level = 0.05, power=0.25) #n needed for small effect, would need n=18 urchins
##pwr.anova.test(k=4, f=0.5, sig.level = 0.05, power=0.25) #n needed for small effect, would need n=4 urchins
##pwr.anova.test(k=4, f=0.8, sig.level = 0.05, power=0.25) #n needed for small effect, would need n=2 urchins